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FORTRAN  PROGRAM  FOR  LINEAR  PREDICTIVE  SPECTRAL 
ANALYSIS  OF  A COMPLEX  UNIVARIATE  PROCESS 


1.  INTRODUCTION 

Spectral  analysis  of  a complex  univariate  process  via  linear  pre- 
dictive and  maximum  entropy  techniques  is  considered  in  reference  1, 
and  Fortran  programs  for  real  data  are  presented  there  in  appendix  J. 

In  this  report,  we  present  a program  for  handling  the  case  of  complex 
data,  yielding  as  an  output  the  auto-spectrum  of  the  process.*  Complex 
data  can  be  encountered,  for  example,  when  a narrowband  real  process  is 
comp lex- demodulated  to  a low  frequency  and  sampled  at  a rate  comparable 
to  the  bandwidth  of  the  process.  When  the  new  center  frequency  is 
zero,  the  process  is  called  the  complex  envelope. 

In  section  2,  an  example  of  the  use  of  the  program  is  presented, 
and  the  changes  that  the  user  must  make  for  his  application  are  pointed 
out.  In  section  3,  the  possibility  of  using  this  program  to  estimate 
the  cross- spectrum  of  two  real  processes  is  investigated  and  found  to 
be  undesirable.  In  section  4,  a limitation  of  the  complex  predictive 
filter  for  complex  waveform  estimation  is  considered,  and  a possible 
generalization  is  indicated  to  alleviate  the  problem. 


*The  theory  and  notation  for  this  case  were  developed  fully  in 
reference  1 and  will  not  be  repeated  here,  for  sake  of  brevity;  the 
reader  is  referred  to  that  earlier  material  for  all  details. 
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2.  USE  OF  PROGRAM  FOR  SPECTRAL  ANALYSIS 

The  program  for  spectral  analysis  of  a complex  process  consists  of 
five  parts:  a main  program  and  four  subroutines,  as  listed  in  appendix  A. 

Input  parameters  to  the  main  program  (listed  in  statement  15)  are 

N,  number  of  complex  data  points, 

PMAX,  maximum  order  of  filter  considered,  and 

J,  size  of  FFT  for  spectral  estimate. 

The  sample  program  generates  a data  example  in  lines  20-33  and  must  be 
replaced  by  the  user  to  fit  his  particular  applications. 

A sample  output  for  N = 100,  PMAX  = 10,  is  presented  below.  It 
indicates  that  PBEST  = 1,  which  agrees  with  the  actual  value  of  p (see 
statement  21  of  the  main  program).  The  fractional  powers  sum  up  to 
0.99999829  instead  of  1;  the  difference  is  a measure  of  whether  the 
spectral  estimate  has  been  adequately  sampled  in  frequency.  (If  the 
error  is  too  large,  J may  be  increased.) 

Since  the  autospectrum  of  a complex  process  is  real,  but  not 
necessarily  even,  it  is  necessary  to  compute  the  spectrum  over  both 
negative  and  positive  frequencies.  Thus  bin  1 corresponds  to  zero  fre- 
quency; bin  J/2  + 1 corresponds  to  +_  Nyquist  frequency,  +^/(2A);  and 
bin  J corresponds  to  frequency  -1/(JA),  where  A is  the  sampling  inter- 
val in  time. 

An  example  of  a spectral  estimate  of  1000  samples  of  a complex 
envelope  of  surface-bottom  forward  scatter  at  a 20°  grazing  angle  at 
frequency  750  Hz  over  a 20  nautical-mile  path  is  presented  in  figure  1, 
where  the  sampling  rate  is  1 Hz.  There  is  observed  to  be  a pair  of 
spectral  peaks  at  +_l/4  Nyquist  frequency,  a strong  very  low  frequency 
component,  and  a rather  symmetric  spectrum  about  zero  frequency.  In  fig- 
ure 2,  the  direct  path  is  employed  instead,  the  sampling  rate  is  0.1  Hz, 
but  1000  samples  are  still  used.  The  center  portion  of  the  spectral 
estimate  reveals  a double  peak  near  zero  frequency  and  a rapid  drop-off 
away  from  this  frequency.  With  these  few  data  points,  resolution  capa- 
bility of  this  quality  is  very  hard  to  achieve  by  any  other  spectral 
analysis  techniques. 
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3.  ESTIMATION  0?  CROSS-SPECTRUM  OF  TWO  REAL  PROCESSES 


Suppose  that  processes  u(t)  and  v(t)  are  real,  zero-mean,  and  sta- 
tionary. If  we  form  the  complex  process 

x(t)  = au  (t)  + Bv(t),  (1) 

where  a and  6 are  complex,  then  the  autocorrelation  of  x(t)  is 

rxx(t)  = x(t)x*(t  - t)  = R*x(-t) 

= ♦ |6|\rv(T)  ♦ a8*Ruv(T)  ♦ a*BRuv(-r),  (2) 

where  the  crosscorrelation  of  u(t)  and  v(t)  is 

ruv(t)  5 u(t)v(t  - x)  . (3) 

The  auto-spectrum  of  x(t)  is  the  nonnegative  real  (noneven)  function 


nxx(f)  = / dr  exp  (-i2irfi)Rxx(T) 


= M R (f)  + |6|\v(f)  ♦ a6*G  i (f)  + a*BG*  (f) , (4) 


w'here  the  cross-spectrum  of  u(t)  and  v(t)  is 


Guv(f)  5 / dT  exP  C-i2TTfT) R„,  fx)  = G*  C-f)  - 


Now  let  us  decompose  the  cross-spectrum  as 


Guv(n  = R(f)  + iI(f)’ 


for  which  (5)  yields 


R(-f)  = R(f),  I(-f)  = - 1(f). 
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Utilizing  (6)  and  (7)  in  (4),  we  obtain 

G (f)  = |ct|2G  (f)  + 1 6 1 2G  (f) 

XX  uu  1 1 vv 

+ 2Re(a6*)R(f)  - 2Im(aB*) I (f) , 

G f-f)  = |a|2G  Cf)  ♦ | B | 2G  (f) 

XX  1 1 uu  1 1 vv 

+ 2Re(aB*)R(f)  + 2Im(aB*) I (f) . (8) 


Solving  (8)  for  R(f)  and  Iff),  the  real  and  imaginary  parts,  respec- 
tively, of  cross-spectrum  G (f),  we  obtain 


Rff) 


I 1 

4 RefaB*) 


[G  (f) 
xx 


G (-f)  - 2 |a | 2G  (f)  - 2 | B | 2G  (f) ] , 
XX  1 uu  1 1 vv 


Iff) 


I 1 

4 ImfaS*) 


Cf) 


G 

xx 


f-f)]. 


(9) 


if  RefaB*)  / 0 and  ImfaB*)  / 0.  Since  a and  B are  arbitrary,  we  choose 
aB*  = (1  - i)/2,  for  which  | a | 2 | B | 2 = | a B* | 2 = 1/2.  For  simplicity,  we 
choose  | a j 2 = j 8 | ^ = \hr2\  then  (9)  becomes 


Rff)  = 


G ff)  + G f-f) 


xx 


XX 


G (f)  + G (f) 
uu  vv 


/2 


= EVENfG  (f))  - 
xx 


G (f) 

uu 


G ff) 
vv 


/2 


Iff) 


G (f)  - G f-f) 
XX  XX 

2 


= ODDfG  (f)  }. 
XX 


(10) 


Now  we  need  only  compute  Rff)  and  Iff)  for  f >_  0,  as  (7)  indicates. 

An  alternative  method  to  (10)  was  presented  in  reference  2,  equa- 
tion (4).  However,  that  method  required  calculating  four  auto-spectra, 
whereas  the  current  method  requires  calculating  only  three  auto-spectra: 
Guuff)  and  Gvv(f)  are  real  and  even,  whereas  I>xxff)  is  real,  but  not 
necessarily  even. 
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The  choices  of  a and  8 in  (1)  are  still  not  unique.  If  we  let 


r1/4ei0,  6 • 2‘1/V4, 


(11) 


then 


-1/2  i(0-4>)  _ 1 - i -,-1/2  -iir/4 


a8*  =2  e 


(12) 


Therefore,  we  must  have  6 - 4>  = -tt/4.  There  are  two  obvious  choices: 
Choice  1 


(13) 


which  is  not  very  symmetric. 
Choice  2 


0 = -tt/8,  <J>  = n/8 


0-l/4  - ir/8  _ 
a=2  e = a - lb 


„ -,-1/4  iir/8 

8 = 2 e = a + lb 


a = 2 cos  (tt/8),  b = 2 sin  (tt/8) 


x(t)  = a[v(t)  + u(t)]  + ib[v(t)  - u(t)], 
which  is  the  preferred  form. 


(14) 


For  known  autocorrelations  or  auto-spectra,  (10)  furnishes  a valid 
way  of  calculating  the  real  and  imaginary  parts  of  the  cross-spectrum 
of  real  processes  uft)  and  v(t).  However,  when  the  auto-spectra  are 
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unknown,  spectral  estimates  must  be  substituted  in  (10).  Although  R(f) 
and  1(f)  can  both  take  on  positive  or  negative  values  in  any  frequency 
range,  there  is  a constraint  on  their  magnitudes.  Namely,  the  magni- 
tude-coherence is  upper-bounded  by  unity.  However,  several  numerical 
examples  using  the  programs  in  appendix  A,  for  cases  where  the  true 
magnitude-coherence  was  near  unity,  yielded  estimated  magnitude-coher- 
ences  greater  than  unity  in  some  frequency  ranges.  This  was  traced  to 
the  fact  that  the  estimate  of  Gxx(f)  can  be  too  small  and/or  the  esti- 
mates of  Guu(f)  and  Gvv(f)  can  be  too  large.  This  type  of  coherence 
estimate  is  intolerable;  hence,  estimation  of  cross-spectra  of  real 
processes  by  means  of  the  auto-spectrum  of  a complex  process  is  dis- 
couraged. The  same  conclusion  is  offered  for  the  method  in  reference  2. 
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4.  A LIMITATION  OF  COMPLEX  PREDICTIVE  FILTER 

The  theory  behind  the  program  presented  in  appendix  A has  been 
given  in  reference  1 and  is  based  on  a linear  predictive  technique. 
Specifically,  given  p past  values  of  complex  process  {x^},  a linear 
one-step  prediction  of  x^  is  attempted  according  to  reference  1,  equa- 
tion 58: 


k ~ ^ anXk-n’ 

n=l 


(15) 


If  we  express  all  the  quantities  in  (15)  in  terms  of  their  real  and 
imaginary  parts  according  to  definitions 

W ‘V  xk  ’ uk  * ivk'  \ * i6n'  (16) 


then  (15)  can  be  expressed  as 

P 


u,  = I (a  u.  - 8 v.  ) , 
k L:  n k-n  n k-n 
n=l 


v = 7 (6  u,  + a v.  ) . 

k L , n k-n  n k-n 
n=l 


(17) 


But  (17)  is  not  as  general  as  the  form  for  prediction  given  by 
(Pn>  vn>  8n>  an  real)  ■ 


u,  = y (p  u,  + v v,  ) , 
k L.  n k-n  n k-n 
n=l 


v,  = y 6 u,  + a v.  ). 
k L , n k-n  n k-n' 
n=l 


(18) 


It  is  apparent  that  the  mean-square  errors  of  predictions  in  (18)  could 
be  made  smaller  than  those  in  (17),  in  general. 


The  complex  estimate  of  xj<  that  can  be  formed  from  (18)  is 

P 

*k  * "k  * i5k  ‘ I/Vk-n  * Vk-n>’ 

n=l 


(19) 
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where 


g„  = + + i(6r,  ' VJ]» 

n l n n n n 


hn  ■ - % * i(6„  * V1' 


(The  form  (19)  is  the  one  alluded  to  in  reference  1,  footnote  to  equa- 
tion (58).)  The  coefficients  (gn)^  and  {hn}^  in  (19)  are  two  com- 
pletely independent  sets  of  complex  constants  that  can  be  chosen.  Since 
the  complex  predictive  filter  in  (15)  is  obviously  a special  case  of 
(19),  it  is  expected  to  have  a more  limited  ability  in  waveform  predic- 
tion than  (19);  however,  (15)  may  suffice  for  spectral  approximation 
purposes.  (For  a real  process,  (19)  reduces  to  (15).) 

Now  suppose  that  a pair  of  real  processes  were  actually  generated 
according  to  p-th  order  autoregressions. 


u,  = y (p  u,  + v v.  ) + d,  , 
k L,  n k-n  n k-n'  k1 
n=l 


v,  = T (6  u,  + a v ) + b,  , 
k , n k-n  n k-n  k 
n=l 

where  all  the  quantities  are  real.  Then  complex  process 


x.  = u,  + iv.  = y (g  x.  + h x*  ) + w,  , 
k k k , lsn  k-n  n k-n'  k 


where  gn  and  hn  are  given  by  (20),  and 


VV  ibk- 


Now  if  an  = -un,  6n  = vn,  for  1 <_  n £ p in  (21),  we  get  gn  = 0,  hn  = pn 
+ ivn,  and  (22)  yields  autoregression 


x,  = I h x*  + w,  . 
k n k-n  k 
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i.  j 


■ 


i 


So  a linear  prediction  according  to  (15)  is  not  expected  to  do  too  well 
on  the  actual  waveform  values  given  by  (24),  no  matter  what  the  statis- 
tics of  (wk)  are.  However,  the  spectral  approximation  available  via 
(15)  of  Gx(f)  for  process  (24)  can  be  very  good  if  p is  chosen  large 
enough  in  (15). 


As  an  example,  let  p = 1 in  (24): 


x,  = hx.*  + w,  , 
k k-1  k 


(25) 


where  h and  wk  are  complex  with 


I h | < 1, 


w,  w*  =6  , w,  w.  = 0. 

k k-n  on  k k-n 


(26a) 

(26b) 


The  excitation  process  described  in  (26b)  is  an  analytic  process,  as 
witnessed  by  the  zero  value  for  the  second  ensemble  average.  We  find 
averages 


k k-n 


and  correlations 


R = x.  x* 
n k k-n 


1 - h 


<R  = x.  x,  - „ 

n k k-n  , i,_i2 


l - |h|‘ 


:-n 

6 , n 

on 

> 0, 

hi". 

n = 0, 

2,  4,  .. 

o. 

n = 1, 

3 , 5 , 

o. 

n = 0, 

2,  4,  .. 

(27) 


(28) 


I h | n~ 1 , n = 1,  3,  5, 


(29) 


Since  (29)  is  not  zero  for  all  n,  process  (xk)  of  (25)  is  not  an  anal- 
ytic process. 


If  we  use  the  facts  (derivable  from  the  definitions  above)  that 

(30) 


R = R*,  (R  = <R  , 
-n  n’  -n  n 
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we  find  that  the  spectra 


G (f)  = A £ R exp  (-i2rfnA) 

n=-oo 


1 ♦ h | 

1 - | h | ^ exp  C-i2irf2A) 


<5(f)  ha  £ f.  exp  (-i2irfnA) 
n=-°° 


2h  cos  (27rfA) 

1 - |h  | 2 exp  (-i2irf2A)  |2 


Generally,  spectrum  G(f)  is  real  and  positive  and  C(f)  is  complex 
and  even  about  f = 0.  For  this  particular  example,  (25),  G(f)  is  also 
even  about  f = 0. 

If  we  attempt  prediction  on  the  process  (25)  according  to  (15)  and 
we  minimize 


we  find 


A | 2 _ | A | 2 

lekl  = lxk  ' X]J 


0,  n i 2 
I h | 2 , n = 2 


■ I A | 2 
min|ek| 


(1  - Ihl2)'1,  p = 1 
1 + | h | 2 , p > 2 


Now  (34)  is  hardly  the  same  result  as  the  actual  autoregression  (25) 
Nevertheless  we  find  spectral  approximation 
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G(f) 


A(1  - M2)'1,  p = 1 
G(f),  p 1 2 


(36) 


from  (34);  that  is,  the  spectral  approximation  is  exact  for  p ^ 2,  de- 
spite the  obvious  error  of  (34)  in  terms  of  waveform  prediction. 

If  instead  we  attempt  prediction  on  process  (25)  according  to  (19) 
and  we  minimize 


(37) 


we  find,  of  course. 


g 


n 


0,  h 


and 


0,n/l 
h,  n = 1 


min 


r?ki2 


(38) 


(39) 


For  | h [ 2 near  1,  the  error  (35)  for  p > 2 is  approximately  twice  as 
great  as  (39) . 

For  the  general  autoregression  in  (22) , it  is  shown  in  appendix  B 

that  for  analytic  white  noise  {w^},  (hn}j  = 0 if  and  only  if  {<Rm Jq  = 

Thus,  given  a complex  data  sequence  {xn}^  of  unknown!  origin,  we  can  de- 
fine 


A 

K 

m 


N 


I 

n=m+ 


x x , 
, n n-m 


0 < m. 


(40) 


Then,  if  |J?m|/^0  <<  1 for  0 m q,  the  autoregressive  model  (19)  with 

{hn}^  = 0 can  be  used  with  some  confidence  for  p < q to  predict  the  ac- 
1 A ~ 

tual  waveform.  But,  even  if  some  f 0,  the  autoregressive  model  (19) 

with  (hn}^  = 0 can  still  be  used  to  estimate  the  spectrum  of  the  process 


15 


TR  5505 


{ X]^ } ; however,  in  order  to  attain  equivalent  spectral  estimates,  it  has 
been  observed  that  more  data  points,  N,.  are  needed  when  some  hn  + 0 

than  when  (hn}j  = 0. 

If  a process  is  generated  according  to  autoregression 


(41) 


where  process  (w]<)  is  not  analytic,  then  «Rm}  are  not  necessarily  zero. 
Thereby  prediction  (15)  will  not  necessarily  give  accurate  predictions, 
although  the  spectral  estimate  can  still  be  adequate;  this  situation  is 
discussed  further  in  appendix  C.  Generation  of  analytic  processes  is 
considered  in  appendix  D,  and  a more  thorough  look  at  the  prediction 
capability  of  (19)  is  considered  in  appendix  E. 
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5.  DISCUSSION 


A program  for  estimating  the  autospectrum  of  a complex  univariate 
process  via  linear  predictive  techniques  has  been  presented.  Although 
it  can  be  used  to  estimate  the  cross-spectrum  of  two  real  processes,  it 
is  not  recommended  because  estimated  values  of  magnitude-coherence 
greater  than  unity  can  result.  Instead,  the  methods  of  multivariate 
techniques  presented  in  reference  3 should  be  employed;  in  fact,  the 
theory  for  complex  multivariate  processes  is  developed  there  and  a 
working  program  given. 

Although  the  program  presented  here  presumes  that  none  of  the  data 
points  are  bad,  it  may  be  readily  generalized  to  include  bad  data 
points.  The  method  and  program  presented  in  reference  1 furnish  the 
necessary  background  for  this  extension. 

Application  of  the  linear  predictive  technique  in  (15)  is  most 
successful  when  the  complex  process  under  investiagtion  is  analytic. 
Otherwise,  the  more  general  prediction  technique  in  (19)  is  worthy  of 
consideration. 


17/18 
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Appendix  A 

FORTRAN  PROGRAM  FOR  SPECTRAL  ANALYSIS 


The  following  program  for  spectral  analysis  of  a complex  process 
consists  of  five  parts:  a main  program  and  four  external  subroutines. 

The  subroutine  BURGCX  computes  the  complex  predictive  filter  coeffi- 
cients, POWERC  computes  the  fractional  power  in  bands  (JA)-1  Hz  wide, 
MKLFFT  effects  a fast  Fourier  transform  (reference  4),  and  QTRCOS  gen- 
erates a table  of  cosine  values  (reference  4) . 
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c SPEciKAL  ESTIMATION  FOR  CO'PLcy  QaTA 
C USl>  I CHANGE  LINE  15  AND  KLPLAcC  LINts  20-33 
C N = i.UMBEr  OF  COMPctX  bATA  POINTS!  IN-j)GER  INPUT 

0 a( i ) , . , . *mim)  = cuiPlex  input  ..atm  Altered  on  ootpui 
C PM A a = MAXIMUM  ORDER  uF  FILTER)  ImTES^I  InPuI 

C J = si/E  OF  FFT  (MUSI  HE  A POER  OF  f TO  USL  . u.FFT ) ) INTEGER  It. PUT 
C PBEST  = BEST  ORDER  OF  FILTER  * INTEGER  OUTPUT 

c au),...,a<PuesT)  = Complex  predict ivt  filter  r 'efficients)  output 

C PROL  = PRODUCT  (1“Al>S(m(P)  )**2)  FOR  P=i  TO  PuEST I uUTpijT 
C RhO(l)»..,,RhO(PMAX)  = COMPLEX  NORMALIZED  CORNEL* T I Oi .SI  OUTPUT 
C XX(1)*...»XX(JJ  = FRACTIONAL  PcWErSI  ot'TPuT 
C COU)»...»CO(D/A+l)  = QUARTER  cOSiuE  TABLE  FOF  FFT  PurPOSES 
E Y IS  A REUUIRLU  COMPLEX  AUXILIARY  ARRAY 

c rr  is  a reouireo  auxiliary  array 

parameter  n=  loo.  pmax=io»  j=  si2.  jai=j/a+i 

INTEGER  PbEST 

COMPLEX  X(N) .Y(n).MPMAx) »RhO(F  1AX) 

DIMENSION  XX(J) .YY(jl.CO( Jei ) 

c complex  input  uata  in  

COMPLEX  Al, 2(14001 

DEFINE  I RAND? I *5**1 5+ ( U-SIr,M < 1 . I*s**lb) ) / 1 > *j4359738367 
DtF INE  RANUrFLOAT ( 1 ) /3A359738387 . 

1=5281 

Al=( .bb, ,85) 

Z(l)=(0.lO.) 
l u 21  L=2»14Q0 
I=IRANb 
Rl=RANb-,5 
1=1RAND 

F<2=HAND-x5 

21  i (L)=Al*2l L-l ) ♦ CMPLX ( R1 » R2 ) 

UO  22  1=1 .N 

22  X<I)=2(IPlA00-N) 

PRINT  1 

1 fokmatuhi.*  input  data:*) 

PRINT  4,  <X(I> . 1=1 *N1 

c evaluate  predictive  filter  coifFicjents 

CALL  BuRGCX ( N . PMAX »X*Y» PBEST. A . PRO, , , RHO ) 

PRINT  9,  X(N) 

9 FoKMAT(/'  KEAN  = (»,L13. 8, *,*,l13. £,*>•) 

R1=R£Al ( Y (N) > 

PRINT  10»  hi 

10  F uRMA  f ( ' STAUUAKO  uEVIatION  = '»r_13.P) 

PRINT  2.  PBEST 

2 FORMAT (/•  FbEST  =*.13) 

IF (PdLST  j£Q . 0 ) 00  10  12 
PRINT  3 

o format ( / * predictive  filter  coefficients  i o«  ).>est:*) 

PRINT  A.  (A(I). 1=1. PBEST) 

A FORMAT (A(£l8<8.Ll5.8iX 

12  PRINT  S.  PROD 
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FORMAT  (/'  PRODUCT  *1-ABS(4(P)  ) ) ='»E13.b) 

PRINT  b 

FuKMAT*/'  NORMAcUlD  CORRELATION  CoEFF IC 1 1 .T s : * > 
PRINT  4,  (RhQ ( I ) » 1=1 »PMAX) 

CALL  QTRCOS ( CO r J ) 

Evaluate  fractional  powers 

call  POwtKC(RbEST»A»PROu.J»xX»rr»Co»SUM) 

PRINT  7 

format*/'  fractional  poaErs:') 
print  b,  ( xx < i i » 1=1 » Ji 
FORMAT(2X*I0E13.fai 
PRINT  11  » SUM 

format*/'  sum  of  fractional  powlRS  z'.eis.G) 

END 


> jlROu  i i nl  i ukulM  is  , PM  A a 1X1  r • PP1'  S r , /■  *PpO!>»  l-h  j) 

Tnis  UUjjRj  rii.t.  Cu-POTtS  TmE  CiU'iPL-  X H-EOICTIV.  ULTLi  COEFFICIENT*. 

N - ,L  ‘.ULk  OF  cOFiPllX  IaTa  POI  TSi  intf  ger  II  I . . T 
p 1A>  - n'AaIMUF  ORoeK  uF  Fl^TLri  IkTEG,  i IuPul 

All)  , XU)  * . . . * A*N)  = LU'iPLi  x Ux.TA  aRRaY  On  IM'.TI  ALT  EKED  ON  OUTPUT 
O.i  OUTPUT,  X * 1 )t  X *2  ),...,  X *PRA>  ) — A * J i Pf"AX  ) , , ( 2 FP  TAX)  t . , , » A *P*’ A>  »F>‘  Ax  ) 
Y*1),Y(2),.,.»Y*N)  = COMPLEX  Ai  X ILI  ARy  APRaYF  SCraTCm  INPUT 
Ok  OUTPUT,  Y ( 1 ), Y * ^ )#,..# Y ( PN«/ ) = A * t ; 1 ) , A ( L i 2 ) , . . . , A IP  AX > PKAX ) 

U,.  luTPUT,  X*N)  = FAiv,  Alip  YU)  = ST„f  OARl  IF'vIaTION  OF  INPUT  DATA 
PnEST  = bL jT  ORDER  OF  PILTrRS  JNTF  (ER  OUTPUT 

A*D ,« *2) , . , . , a<polST)  = complex  r k£d i c t i /l  filter  Coefficient  array  r 

A * 1 J PuLST  ) , A * 2 ! PbtST  n ( Pur  ST  F P3£t,T ) F O'  1PUT 

Pkui  = PRu'  UCT  ( l”MbS  ( « ( P J p ,,E  ST  ) ) **c  > pcR  P=1  TO  r uES  T J OUTPUT 
RmOU),...,RHO*PF,ax)  = COMPLEX  NO r ALj/ED  COR^lla t IC; IS t OUTPUT 
CuMi-llA  X ( k j , Y (N)  , A *PwAX ) , i-OiCj  * pMAX  ) jS  RECUIPE'i  I H r.AlN  pROfoReM 
INTEGER  PMA/»PueST»P 
LvjOUlE  . KlCISIOi  SmK»SA1»SI: 

C Ui*'PlE  a SI » b » T 

COMPLEX  X*1),Y(1),«(I),  r.HO  ( 1 ) 

IF(P'TAX.GT.3.*bo/RT  *(.) ) pRInj  2»  P.Ma>»N 

FvikMaT*/'  Pi-ax  = * » 14  » • IS  TrO  LaRGe  FOR  <UMP->  OF  LATA  POINTS  II  =' 
4.1b) 

COvFofE  I'LAIj 

jl- ( ) • , C • ) 

V J 1 I = ) i N 
L I — S 1 MX ( I ) 

S 1=CMPlXFHEAL(S1)/N,AIi»MG*SJ  )/N) 

SnuTkaCT  L An » aNu  SCmLE  Tb  UNjT  vARI/NCE 

bc=0. 

iJ  J 1=1 *N 


> ( I >=X  * I )-Si 

S«-  = S2+rEAL*X*1)  ) **£  + A IMaG * X ( I ) ) * *2 
S<-=SOHT  (Si./  (N-l  . ) ) 
i=l# /it 
L J 5 1 = 1 »N 

a i * ) =C..,r  LX*r-EAL*X*  1)  )1»B,  A IMAGO  ( I)  )*H) 
I ( 1 ) =X ( ( } 
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r 


i 

> 

! 


i 

- 


I 


L 


i Jr.'Ji S 1 Cm  . 

rrt 

r n \JuJ C — 1 • 

» ioMNsr» 
i t>LST —0 
t hyL  — 1 ■ 

0 i ;i'*  1 

1 c*i_cuy.*TL  .-RObb-GAini  f_u.  15b 

1.IIIC1;  .U"1 

ShI  - / #U#' 
bdi:>.L>t 

l-F*1 

i o 7 I =i_  * •« 

1 1=KEAL< XI I ) 1 
l.-=Ali».Ar-tA(i)  ) 

r o=«<Lah  if  1 1-1  > ) 

1a  — Alf-'AvilYli  — 1)) 
jAKiS'W\  ♦■Tl  * T.S4  ( *:*  T o 

i*Sm1  ♦ 1 11*  t r l*  i ^ 

7 bo=Sb+Tl**2*T2'**2  + T;5**24-T4**2 

l=^.*SAk/Sb 
C=2.*SAI/StJ 
G=C  4PUX(B#C) 

L=l.-d*B*C*C 

PK00UC=PK0UUt*d 

o CALCULATE  FILTES  CuEFF’ICIEl i TS * EUS.  STORE  IN  X (1 J . . . , , X < P) 

X(P)=G 

IK(P.Ey,,lJ  GO  Tu  a 
L=P/2 

LO  9 I=lfL 

T=A(l)-G*CONdG(X(P-ll) 

A (P-I ) =X  <P"I > -G*COhJG( A ( I ) ) 

9 X(I)=T 

C CALlULATE  nOKMALI2ED  CORRELATION  COEFFICIENT!  FV.  199 
o T=A(P) 

IF(P.EG.l)  GO  TO  19 
L=P-1 

LO  lb  I r 1 » L 

lb  1=T>A(l)*KhO(P-I) 

19  AhO(P)=T 

C CALCULATE  aKAIkE’S  INFORMATION  CRITERjGNI  EoS,  ltb«202 
PELERHzB*bNGL(S9)/(2.*(i.-P)  ) 

AlC=LUG(RLLL8Ki+<*.*FLUAT  <P)/(N-PJk 
IF  IAIC.GE.AICMIn)  00  TO  10 
AlCMlNzAlC 
PdESTzP 
FKOUzPNGUUC 
LO  11  I=1»P 
11  >.(I)=X(I1 

10  IF IP.Eo ,PMA* ) oU  TO  lb 

L UPumIE  FOK^ARO  ANU  UACKftAKo  SEQUENCES,  EG.  lo3 
L=P+1 

to  12  I=N»L*-1 

T=a ( I )»G*Y ( 1-1 l 

Y(l)sY(Ivl)-CONjG(bj*X(l) 
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12  A(I)=T 
YiP)=G 
00  TO  b 

lb  Y(PMAX)=6 

IF  (PriEST pEG.PMAXl  oO  TO  4 

i c compute  extrapolated  normalize^  cci-RELATin:, 

c coefficients  from  pbest+i  to  fvaxi  Eu.  16s 

L=PbES  T ♦ ! 

10  17  P=L»PMAX 
AiP)=(0, »0. ) 

T=(0.»0.1 
LO  Id  1=1 »PbEST 
IB  1=T+A( I)*KHO(P-I) 

17  prto(P)=T 

4 a (N) =S1 

Y(N)=CMPLX(S2»0.) 

aeT  URN 

C.14O 


SuuKOU  I Iul.  pO*l.i<C  IrLEST  * A »l'i.  OU» X/ » YY  » Cl  . SU  ) 

TniS  bUdRui  TINE  Co-IPOIlS  ThE  FuACT  10N,,L  Pu.it  S ii.  HAi.,,S  1 / ( J*UI  LTa) 

P.JESI  = Ut-'oT  OpDEn  oF  FILTER)  INTEGER  INPUT 

m (1 ) , . . . r a ( PbLbT  ) = Col  PlEa  HLTEb  CO(  FPICli  !•  T i.mAY)  IIVuT 

PkUL  = PkuuUCT  ( 1-AoS  (h  ( P)  ) **2  ) FOt  P=j  TO  F'EbT)  INPUT 

U = SIZE  0(  FFT)  I lTEoLi<  IrlPUl 

Xa  = AUXILIARY  ARnaY  olJ  INPUT 

XX(  1 )»...» XX(u)  - FRACTIONAL  Kjft'Ef  E 0;.  OUTPUT 

YY  = AUXILIARY  ARRAY  I SCRATCh  [NPuT 

Colli»...fC0tU/4-H)  = OUARTER  c°Si,.E  t,U)Ll  Fo.v  mT»  IiP"T 
OlMLnblon  /Alu)  r Yr  ( JJ  ,cu(  U/4+1  ) 1'..  Rfi.UlREO  1 1 aIP  pi<Oir<AM 

CO**1.  w^-X  A | 'MAX  ) 1 j Ki_  OUlntD  J'  M/-  I ’)  | F OGn A ’ * » itiiF  ) -AX*  OE  • i • 1 1*  S | 

iulEGLk  PbLjT 

uI-iENSl:;N  xX(  1)  , Y Y C l»»Co(l> 

COMPLEX  All) 

F=PROO/.J 

)Xll)=l. 

Y Y 1 1 ) = C . 

IF  IPuL-jT.tu.C  ) jO  10  4 
Lu  1 1=1 .PLLST 
> a ( 1*1 ) =tRlaL  l h l 11) 

Y 1 l I + l ) r-AlI-  AO  l A ( 1 ) ) 
l_=PbL5T  + 2 

L 0 2 I=L*U 
XX  ( I )=(.', 

YY (I)=C. 

0=1.442 ?*uOo(J)+.b  U LoGPl.n 

C.«LL  i-'im.F  F T ( Xx  • Y Y » LC  * L » - 1 ) 

SUi-.=  T. 
iii  i 1 = 1 # J 

/Ail) =f / < xx ( 1 ) * *2 ♦ir(I)» ♦ <? ) 

S J'’:  = o0l  4 XX  ( I ) 

1 UJMI. 

L 4^ 
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bvo'^OUT  INL  i'.IS  LFr  T ( A » Y * CO  # M » I SN) 
l i>i£(JSlCN  X(1 ) *Y (1) »CC<1>  »L(X2J 

l o-wIvALtNtt  (Lli.Ltl)).  (Ll),L(i;))»(llO»L(j))MLt),Lll4))t(Uii  ("i)). 
XlL7»Llb)X*(lb>L(71),»Lb,LtP)>»lL4.L(9)),(L3,Llln)>.<L**L<lX)>, 
2li_i*LtX2)  * 

l .j4=U/4 
i.u^HUlvH+l 
I L4F2  = l»j4Pl4l 
>g2P2=.'Jg4+M4P2 
l g b LUrX  # i*i 
L>iA=2**  (fc-LC) 

LIX=2*L  'X 

lbCL=.i/LiX 

tU  <1  Li^XfLXX 

lAnG=(L'-rlJ*ib(.L+l 

IF  11ARG, LL , hOLP l ) gC  To  4 

C=“CC  (i-*02P2-lAr<b) 

b=lb)i*CC  « IAHG-Imu4) 

GU  TO  o 
*4  (.=CC(ImFG) 

j=iS,:*Cf  «NL4R2-1aLu) 
o lg  b LI -LI X t U» LiX 
ul=Ll-LlX+LX 
v2=J1+L''X 
1 1=X ( Ji )»X(02) 

T«i=Y(Jl)*Y  (02> 

AibX)=X(Ui>+XTj  c.) 

Y (wX)=Y (jX)+r (Jt) 

Xib2)=C*TX-S»T2 
Y (G2)=C*T2+S*T1 
o LghTINUF 
LG  4C  J-X ’ X2 
HO)=l 

IF  ( J-*’ ) JX*3X*4C 
Ox  L l O ) — i.*  * T F'4  X” b ) 

4C  Cgi.TINl/p 

v-4— 1 

L v bf  J1?X*LX 
iv  t>"  *j2?OX»L2»l1 
t u oC  g3»g2#L3»i_2 
i 0 br  v'‘=03*L4m.0 
Lo  b'  jb; J4 f Lb  * l4 
v v bC>  J6?bb»Lfa»Lb 
LG  br  j7SGb*L7»Lb 
vb  b0  JSsOTfLbu.? 
uU  o'  v9=vo»L9»Lb 
iu  bT  u1Q=J9»LX3»L9 
i »g  b?  JlX=JX0#LXX»Ljn 
Lg  bC  jfcJXX t LX2 • LX X 
IFlJ(-JKX  bX*bX,b2 
51  Fsa(JN) 

> (Gll)=A  ( JKJ 

> l vR ) SK 
F l-t  (JtO 
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6C 


Y < JIJ)=Y  ( JK) 

Y ( JN)=H 
JhSJlM+1 
CONTINUE 
KE1  Ulvf* 

L.Jb 


^uoKOUTlNt  (jTKCoSUm) 
L li-itiiSiUN  t (1) 

, >+1=,j/<4  + 1 
bLi-=b.^tjiit3b3C7/,J 
10  1 I=1»N41 

t ( n=cos<  < i-i  )'*scu 

huT  Jr<f  i 

t.lMU 


) 
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Appendix  B 

PROPERTY  OF  AUTOREGRESSIVE  MODEL 


The  process  of  interest  here  is  given  by  the  autoregressive  model 


(22) 


x,  = y (g  x.  + h x*  ) + w,  , 
k L i V6n  k-n  n k-n7  k’ 
n=l 


(B-l) 


where  excitation  (w^)  is  analytic  white  noise.  That  is, 


w,  w*  = 6 , w,  w,  = 0 , 

k k-m  om  k k-m 


(B-2) 


It  then  follows  easily  that 


w.  x*  = 6 , w.  x.  = 0 , m > 0 . 

k k-m  om  k k-m  ’ — 


(B-3) 


Use  of  (B-2)  and  (B-3)  then  yields  correlations 

P 

R = x x,  = y (g«  + h R*  ) , 

m k k-m  &n  m-n  n m-n 

n=l 


m > 0, 


( B— 4) 


R = x,  x* 
m k k-m 


y (g  R + h R*  ) + 5 , 

&n  m-n  n m-n  om 

n=l 


m > 0. 


(B-S) 


For  given  coefficients  (gn)^  and  (hn}^,  (B-4)  and  (B-5)  constitute  si- 
multaneous equations  in  the  unknown  correlations. 

Now  let  us  suppose  that 

h^  = 0,  1 _<  n _<  p . (B-6) 

Then,  from  (B-4),  the  first  p + 1 equations  yield 


B-l 


; ; 


I ' 


i 

-v 


I i 


■ 
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®o  - glRl  - • 


. - g f = 0 
P P 


% ■ giVi 


g K = 0. 
p 0 


Therefore 


I?  = 0,  0 < m < p, 
m — - 


and  from  (B-4)  and  (B-6)  it  follows  that  is  zero  for  all  m. 
Conversely,  assume  that 


R = 0,  0 < m < p. 
m — — 


Then,  from  (B-4),  the  first  p equations  yield 


h R + h R + ...  + h R = 0 

11  2 2 p p 


hiR2-P  * h:R3-p  * * Vi  ’ °' 


Therefore 


h = 0,  1 < n < p, 
n _ _ 


(B-7) 


(B-8) 


(B-9) 


(B- 10) 


(B-ll) 


and  from  (B-4)  and  (B-ll)  it  follows  that  (Rm  is  zero  for  all  m. 

Thus  we  have  shown  that  {hn)P  = 0 if  and  only  if  = 0 in  the 

autoregression  (B-l)  with  analytic  white  noise  excitation. 


B-2 


. 1 
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Appendix  C 

NONANALYTIC  WHITE  NOISE  EXCITATION 


Suppose  a process  is  generated  according  to  autoregression 


xk  ' !,Vu  * V 

n=l 


(C-l) 


where  excitation  (w^)  satisfies 


w w*  = B6  w.  w,  = 66 

k k-m  o rri  k k-m  om 


6 is  complex  and  nonzero;  therefore  {w^}  is  not  an  analytic  process, 
although  it  is  white. 

Then  {Rm}  need  not  equal  zero,  even  for  autoregression  (C-l).  For 
example,  for 


P = If  gx  = g»  I g I < 1,  g complex. 


we  find  correlations 


r = — - — 2 > Rm  = gV» m > i>  r m = r:  , 

O ^ I g | m 0 _ ~m  m 


<R  = 6 , <R  = gmK  , m > 1,  f?  = R . 

o , 2 ' m b o’  — -m  m 

1 - g 


The  corresponding  spectra  are 


G(f)  = 


|l  - ge 


-i27TfA,2  ’ 


C(f)  = 


1 - 2g  cos  (2irfA)  + g 


(C-3) 


(C-4) 


(C-5) 


(C-6) 


(C-7) 


C-l 
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• - — 


I 


f 


Since  the  {Rm}  are  not  zero,  appendix  B shows  that  the  coeffi- 
cients {hn}  would  not  be  zero  in  model  (B-l)  with  an  analytic  excita- 
tion. 


Optimum  prediction  on  process  (C-l)  using  (15)  gives 


(C-8) 

with  a minimum  mean  square  error  equal  to  B,  and  the  spectral  estimate 
is  identically  (C-6) . Thus,  for  this  example,  the  nonanalyticity  of 
the  excitation  is  no  problem. 

For  the  more  general  model  of  (C-l)  with  p > 1,  it  can  be  shown 
that  all  { Rm } are  independent  of  the  value  of  6.  Then,  although  (15) 
may  not  be  too  accurate  for  waveform  prediction,  it  can  still  be  used 
for  spectral  estimation  purposes. 

Optimum  prediction  on  process  (C-l)  using  (19)  gives 

g^  = g,  other  coefficients  = 0,  (C-9) 

with  a minimum  mean  square  error  equal  to  B.  This  yields  the  same  re- 
sult as  above. 
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Appendix  0 


A METHOD  OF  GENERATING  ANALYTIC  PROCESSES 


Suppose  linear  filter  H(f)  is  excited  by  complex  input  x(t), 
yielding  output  y(t) . Then  correlation 


R (t)  H y (t)y*  (t  - t)  = / df  exp  (i2irfr)Gx(f)  |H(f)  | , (D-l) 


and  spectrum 


G (f)  H / dt  exp  (-i2irfi)R  (x)  = G (f) |H(f) | . CD-2) 

y y *■ 


Also  correlation 


<Rv(t)  = y(t)y(t  - t) 


= / df  exp  (i27ifi)Cx(f)H(f)  H(-f), 


CD—  3) 


and  spectrum 

C (f)  = / di  exp  C-i2rrfT)R  (t)  = C (f)  H(f)  H(-f)  . (D-4) 

y y * 

Complex  process  y(t)  is  defined  as  being  analytic  if  CD- 3)  is  zero 
for  all  x.  Suppose  that  filter 


H (f ) = 0 for  f < 0. 


(D-5) 


Denote  the  output  of  filter  (D-5)  by  y+(t).  Then  (D-4)  shows  that 
<?v(f)  and  P v(t)  are  both  identically  zero  for  all  argument  values. 
Therefore  single-sided  waveform  y+(t)  is  an  analytic  process. 

Let  complex  envelope 

£(t)  = y+(t)  exp  (-i2irfot).  ( 


D-l 


I 


i 
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R (t)  2 y(t)y*(t  - x) 


R (-0 

y+ 


for  any  fQ.  Thus  the  complex  envelope  of  any  stationary  process  is  an 
analytic  process. 

On  the  other  hand,  for  the  two  real  processes  u(t)  and  v(t),  no 
linear  combination,  au(t)  + Bv(t),  where  a and  & are  complex,  ever 
yields  an  analytic  process  unless  Ruu(t)  » Rvv(T)»  anc*  RuvCT)  satisfy 
very  special  restrictions.  Thus  the  process  constructed  in  (1)  was  not 
analytic  and  could  not  have  been  expected  to  yield  good  prediction  cap- 
abi lity  via  (15) . 
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Appendix  E 

CAPABILITY  OF  A MORE  GENERAL  PREDICTION  MODEL 


If  p is  infinite  in  (19),  we  have  prediction 


x,  = 7 (g  x,  + h x*  ) . 

k L , n k-n  n k-n 
n=l 


Minimization  of  |e*jJ2  - Ix^  * x]<  | 2 yields  the  simultaneous  equations 


7 (gR  + h f * )=R,l<m, 

L , 6n  m-n  n m-n'  m’  — ’ 

n=l 


y (g  <R  + h R*  ) = <R  , 1 < m. 
n m-n  n m-n  m — 
n=l 


It  can  then  be  shown  that  is  a white  process  with 


I ~ | 2 

mm  e. 


= R - I (g  R*  + h «*). 
o u , n n n n 
n=  1 


Also  it  can  be  shown  that 


t.  eT,  = 0 for  m / 0, 
k k-m  ’ 


e,  = K - X Cg  * h R 
k o L . n n n n' 
n=l 


However,  is  not  an  analytic  process  since  (E-5)  is  not  zero. 
The  simplest  example  to  demonstrate  this  is 

R=R  6 , K = E 6 , (E 

n o on’  n o on’ 

for  which  (E-3)  and  (E-5)  yield  R0  and  <R0,  respectively. 
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The  spectral  relations  for  (E-l)  take  the  forms 


G*j»(f)  = Gx(f)|A(f)|2  + Gx(-f)|B(f)|2 


+ Gx(f)  A (f) B*  (f)  + [Cx(f)A(f)B*(f)]* 


and 


Vf) 


«x(f)A(f)A(-f)  + C*(f)B(f)B(-f) 


+ Gx(f)A(f)B(-f)  + Gx(-f)A(-f)B(f), 


where 


and 


A(f)  = J £ exp  (-i2irfnA) 
n=l  n 


B(f)  = ][  h exp  (-i2rfnA) 

n=l  n 


(E-7) 


(E-8) 


(E-9) 


are  considered  known  after  solution  of  (E-2)  for  coefficients  (gn)  and 
{hn}.  Equations  (E-7)  and  (E-8)  can  be  solved  for  Gx(f)  and  Cx(f) : 
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A (-  f ) , 

4 

c,  = 

X 

Gx(f), 

G 

X 

" ('x("f)’ 
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(E- 10) 


(E-l  1) 


E-  < 
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This  requires  the  inverse  of  a 4 x 4 complex  matrix  at  each  value  of 
frequency  f. 

Thus  an  alternative  spectral  estimation  technique  is  available 
from  the  more  general  prediction  model  in  (F.-l).  Whether  it  is  worth- 
while in  terms  of  stability  and  resolution  is  unknown,  as  it  has  not 
been  pursued. 
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